# ==============================================================================
# wolves_map_trend.R
# author: Anselm Hager / Bernhard Clemm
# ==============================================================================

dir <- dirname(dirname(rstudioapi::getSourceEditorContext()$path))
source(paste0(dir, "/code/setup-packages.R"))

# MAP ==========================================================================

# Load wolf data: aggregated number of attacks per municipality
wolves_by_municip <- read.csv(paste0(dir, "/data/wolves/wolves_by_municipality.csv")) %>%
  mutate(ags = ifelse(nchar(ags) == 7, paste0("0", ags), ags))

# Load 2021 shape file ####
shape <- read_sf(paste0(dir, "/data/crosswalks/shapes_2021/vg250gem.shp")) %>%
  rename_all(tolower) %>%
  mutate(bundesland = stringr::str_extract(ags, "^.{2}")) %>%
  rename(gem = gen) %>%
  select(ags, gem, bundesland)

# Merge wolf data with municipality shapes  ####
shape <- merge(
  shape, wolves_by_municip %>%
    select(ags, attacks),
  by = "ags", all.x = T
)

## kick out states for which we don't have data: MV ("13"), RLP ("07")
missing_states <- c("07", "13")
shape <- shape %>%
  mutate(state = substr(ags, start = 1, stop = 2)) %>%
  mutate(attacks = ifelse(is.na(attacks), 0, attacks)) %>%
  mutate(attacks = replace(attacks, state %in% missing_states, NA))

# Categorize number of attacks ####
shape <- shape %>%
  mutate(attacks_fac = factor(case_when(
    attacks == 0 ~ "0",
    attacks %in% 1:5 ~ "1-5",
    attacks %in% 6:10 ~ "6-10",
    attacks > 10 ~ ">10"
  ),
  ordered = T,
  levels = c("0", "1-5", "6-10", ">10")
  ))

# Plot  ####
pal <- brewer.pal(9, "YlOrRd")

map_plot <- ggplot() +
  geom_sf(data = shape, aes(fill = attacks_fac), lwd = 0.05) +
  scale_fill_discrete(
    na.value = "grey70", type = c("white", pal[4], pal[6], pal[8])
  ) +
  theme_minimal() +
  guides(fill = guide_legend(title = "Attacks")) +
  theme(
    legend.position = "bottom",
    legend.key = element_rect(color = "grey30"),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    plot.margin = margin(r = -0.3, l = -0.3, t = 0.3, unit = "cm")
  )

# TREND ========================================================================

# Read in data ####
## Since we do not need municipalities for the overall trend,
## in this data set we have integrated yearly aggregates for
## (1) cases from Schleswig-Holstein after 2020
## (2) cases from Mecklenburg-Vorpommern
## (3) cases from Rheinland-Pfalz

wolves_by_year <- read.csv(paste0(dir, "/data/wolves/wolves_by_year.csv"))

# Plot ####
trend_plot <- wolves_by_year %>%
  ggplot(aes(x = year, y = value, group = name, fill = name)) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_y_continuous(breaks = c(0, 1000, 2000, 3000, 4000)) +
  scale_fill_manual(
    values = c("grey70", "grey40"),
    labels = c("Attacks", "Killed animals")
  ) +
  theme_light() +
  theme(
    legend.position = "bottom", # c(0.3, 0.8),
    legend.title = element_blank(),
    axis.title = element_blank(),
    legend.background = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.text = element_text(size = 14),
    axis.text = element_text(size = 12),
    plot.margin = margin(r = 0, l = 0, t = 1, unit = "cm")
  )

# COMBINE FOR FIGURE ===========================================================

ggarrange(trend_plot, map_plot,
  widths = c(1, 1.5)
)
